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Abstract. The structural and dynamic properties of silica melts under high 
pressure are studied using molecular dynamics (MD) computer simulation. The 
interactions between the ions are modeled by a pairwise-additive potential, the 
so-called CHIK potential, that has been recently proposed by Carre et al. The 
experimental equation of state is well-reproduced by the CHIK model. With 
increasing pressure (density), the structure changes from a tetrahedral network 
to a network containing a high number of five- and six-fold Si-O coordination. In 
the partial static structure factors, this change of the structure with increasing 
density is reflected by a shift of the first sharp diffraction peak towards higher 
wavenumbers q, eventually merging with the main peak at densities around 
4.2g/cm^. The self-diffusion constants as a function of pressure show the 
experimentally-known maximum, occurring around a pressure of about 20 GPa. 



PACS numbers: 61.20. Ja,64.70.ph,66.30.H-,62.50.-p 

1. Introduction 

Under ambient conditions, amorphous silica (Si02) forms a disordered tetrahedral 
network. Tetrahedral units with silicon atoms in the centre are connected to each other 
via the oxygen atoms at the four corners of each tetrahedron. Under compression, the 
network structure of Si02 may change significantly. From experiments on silica glasses 
[TJ[21[31[ll[Sl[ni[7], it is known that the SiO coordination changes for pressures of the 
order of 15-20 GPa. Then, the network structure is dominated by five- and six-fold 
coordinated SiO units. Albeit of great interest for geoscience, experimental studies 
of pure silica melts under (high) pressure are very rare [6]. This is due to the high 
temperatures that are required to study silica in its liquid state. Hence, computer 
simulations are an important tool to elucidate structural and dynamic properties of 
silica melts at high pressure. 

In many of the recent molecular dynamics (MD) computer simulations, the pair 
potential proposed by van Beest, Kramer and van Santen, the so-called BKS potential 
[8], has been used to model the interactions between the Si and O ions. Tse et al 
[9] showed that in BKS silica the oxygen coordination around silicon atoms changes 
from 4 to 5 at intermediate pressures (10-15 GPa) and reaches 6 at high pressures. 
Subsequent studies of BKS sihca [101 El [12] have indicated that transport processes 
first become faster with increasing pressure, in agreement with experimental studies 
on silicate melts [13l [M] [15l [16] . In Ref. [10] , this was interpreted as a signature of a 
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strong-to-fragile transition, according to Angell's classification of glassforming liquids 
[Tfl ITS] . Note that the anomalous transport in silica under pressure was already seen 
in a pioneering MD simulation by Angell et al [12] using a different model potential. 
Further studies on BKS silica [501 HD 123] addressed the question to what extent 
the properties of silica are typical for a strong glassformer, sharing its properties with 
other network-forming liquids such as water. In Ref. [20j . the possibility of a phase 
separation between a low and high density liquid phase was discussed and a rough 
estimate of the free energy as a function of volume and temperature was used to 
locate such a transition for BKS silica, finding a critical point at a density of about 
2.9g/cm^ and a temperature of about 2000 K. 

But recent works have called in question that simulations with simple pair 
potentials can give reliable results for pressure-induced phenomena in silica. Trave et 
al [25] claimed that empirical force fields such as the BKS model "are not sufficiently 
accurate to reproduce pressure-induced phenomena such as changes in coordination 
and increasingly strained geometries" . These authors used "first-principles" Car- 
Parrinello molecular dynamics (CPMD) simulations to investigate liquid silica for 
pressures up to about 30 GPa at a constant temperature of T = 3500 K. In a CPMD, 
the electronic degrees of freedom are taken into account in the framework of a density 
functional theory |24j and thus it is not relying on the use of effective interatomic 
potentials. Whereas the CPMD calculation for silica shows good agreement with 
respect to the experimental equation of state, the BKS model underestimates 
the experimental specific volume systematically by about 15%. Based on CPMD 
calculations, Tangney and Scandolo (TS) [26] have developed a new effective potential 
that nicely reproduces the experimental equation of state of amorphous silica as well 
as some of the properties of crystalline phases The TS model is a fluctuating 

dipole moment potential, accounting for the polarizability of the oxygen atoms. A 
disadvantage of the TS model is its low efficiency. Due to the expensive self-consistent 
computations of the dipole moments on the oxygen atoms in every time step, up to 
two orders of magnitude more computer time is needed in MD simulations with the 
TS model than for comparable simulations with a simple pair potential such as the 
BKS model [57]. 

The conjecture of Trave et al about the validity of simple effective interaction 
models for silica has been recently challenged by Carre et al [28] who developed a 
new pair potential for silica, the so-called CHIK potential. The CHIK potential was 
also parametrized from CPMD calculations and was shown to be superior to the BKS 
model with respect to various properties of amorphous silica (in particular the density 
at low pressures) [28]. In this work, we show that it also reproduces very accurately 
the experimental equation of state, various structural properties and the anomalous 
diffusion dynamics of silica under pressure. This proves the transferability of the 
CHIK potential to the description of amorphous silica under high pressure. Thus, the 
CHIK model can be used to readdress the questions about the possibility of a liquid- 
liquid phase separation and the mechanisms of anomalous transport phenomena at 
high pressures. 

2. The model potential and details of the simulation 

At first glance, it might be surprising that the directional covalent bonding in silica 
can be described by an interaction model that depends only on the magnitude of the 
distance between atom pairs. In fact, in a one-component system such as silicon one 
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Figure 1. Pressure as a function of density for several temperatures, as indicated. 
The experimental data are taken from Ref. |33) . 



has to include at least three-body interactions to reproduce a tetrahedral network 
structure. However, in a binary system such as silica there is the competition 
between three different pair potentials. By tuning the parameters for the repulsive 
00 interactions and the attractive SiO interactions appropriately, one can generate a 
network structure with Si-0 tetrahedra as the preferred structural units. 

Apart from a rather unreliable description of the equation state, the widely 
used BKS potential [5] has been proven successful to describe various static and 
dynamic properties of amorphous silica quite reliably (see, e.g., Ref. [301 131] )• 
The paramctrization of the BKS model was based on Hartree-Fock calculations of a 
single Si04 tetrahedron, charge-saturated by four hydrogen atoms [5]. In contrast 
to that, the development of the CHIK potential was based on a CPMD simulation 
of a bulk system of 114 particles. A structural fitting procedure was used for the 
paramctrization. The idea was to match the partial pair correlation functions, as 
obtained with the new effective potential, with those obtained from CPMD. As the 
BKS potential, the CHIK model consists of a long-ranged Coulomb part and a short- 
ranged Buckingham potential, 

Uap{r) = — 1- Aap exp {-Bapr) (1) 

where r is the distance between an ion of type a and an ion of type (3 (a, [3 — Si, O), 
e is the elementary charge, and gsi — 1.910418 and qo = —qsi/2 are effective partial 
charges. The values of the parameters {Aa/s, Ba/SjCap} are ^sisi = 3150.462646 eV, 
Bsisi = 2.851451 A-i, CsiSi = 626.751953eV/A6, Agio ^ 27029.419922eV, 
Bsio = 5.158606A-1, Csio = 148.099091 eV/A^, Aqo = 659.595398eV, Bqo = 
2.590066A-\ and Coo = 26.836679eV/A6. More details on the potential can be 
found in Ref. [28]. 
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Figure 2. Partial pair correlation function for the Si-O correlations, (;sio('')i for 
different densities at the temperature T = 2750 K. The arrow indicates the cut-off 
''cut = 2.35 g/cm'^, used to identify SiO bonds. The inset shows an enlargement 
of the data. The two vertical dashed lines at 1.618 A and 1.647 A indicate the 
location of the maxima for the densities p = 2.3g/cm^ and p = 3.9g/cm'^, 
respectively. 



MD simulations have been performed for systems of iV = 1152 particles at the 
densities p = 2.3, 2.4, 2.5, 2.7, 2.9, 3.0, 3.1, 3.3, 3.5, 3.6, 3.7, 3.9, and 4.2g/cm3 in the 
temperature range 6100 K> T > 2100 K. The equations of motion were integrated by 
the velocity form of the Verlet algorithm [32] using a time step of 1.6 fs. The long- 
range Coulomb forces and energies were computed by the Ewald summation technique 
[32]. At each state, the system was equilibrated in the NVT ensemble during 0.2 ns 
to 25 ns. Temperature was kept constant by coupling the system periodically to a 
stochastic heat bath. After the equilibration, structural and dynamic quantities were 
obtained from production runs in the microcanonical ensemble. 

Figure [1] displays the equation of state, p{p), along different isotherms. For 
pressures up to about 10 GPa, the pressure is almost independent of temperature in 
the range 4300 K> T > 2100 K. In this pressure range, the simulation is in very good 
agreement with experiment [33] . This demonstrates the transferability of the CHIK 
model: this model has been parametrized from a single CPMD run at T = 3600 K and 
p = 2.2g/cm'^, without considering states at high pressure. Thus, at least as far as 
the equation of state is concerned, the CHIK potential provides a reliable description 
and so it appears to be appropriate for the investigation of silica under high pressure. 

3. Results 

To analyze the change of the network structure of amorphous silica under compression, 
we consider the temperature dependence of the Si-O as well as the 0-Si coordination. 
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Figure 3. Temperature dependence of different SiO and OSi coordination 
numbers z. a) p = 2.5g/cm^, b) p = 2.7g/cm^, c) p = 3.1g/cm^, d) 
p = S.Sg/cm"^, e) p = 3.7g/cm^, f) p = 3.9g/cm''. 
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Figure 4. Pressure dependence of the Si-O coordination at T = 3580 K, 2230 K, 
and 2100 K. 



To this end, coordination numbers zgio and zqsi are defined as the number of O (Si) 
atoms surrounding a Si (O) atom within a distance r < Tcut- For the cut-off radius we 
use the value rcut — 2.35 A that corresponds to the location of the first minimum of the 
partial pair correlation function for the Si-O correlations, gsio(''') [HI- As Fig.[2]shows, 
the location of this minimum is indeed essentially independent of density. However, 
one can infer from the inset of Fig. [2] that the length of the SiO bond changes slightly 
from about 1.62A at 2.3g/cm^ to about l.fisA at 3.9g/cm'^. This finding is in very 
good agreement with results from X-ray diffraction experiments 4J. Moreover, the 
periodicity of the peaks in 5sio('') changes with density. As we shall see below, this 
is due to the loss of tetrahedral order at high densities, reflected, e.g., in the Fourier 
transform of (7sio('')j the partial structure factor S'sio(9): by a shift of the first sharp 
diffraction peak to higher wavenumbers q. 

Figure [3] shows the probability distributions P{z) at different densities for all the 
relevant coordination numbers zgio and zqsi as a function of temperature. For the 
three densities p = 2.5, 2.7, and 3.1 g/cm^, most of the Si atoms exhibit the tetrahedral 
four-fold coordination by O atoms and most of the O atoms serve as bridging atoms 
between the tetrahedra (zosi = 2). In this case, the local defect structures that 
correspond to zsio = 5,6 as well as zosi — 3 become less frequent with decreasing 
temperature such that one expects a very low probability of these structures at low 
temperatures. A different behavior is seen at the density p = 3.5 g/cm^. Now, zsio = 5 
has a higher probability than the tetrahedral coordination 2:310 = 4 and one observes a 
relatively high probability of zqsi = 3 and a significant number of six- fold coordinated 
Si atoms. One can hardly extrapolate the curves for 3.5g/cm'^ to lower temperatures 
since at this density all the coordination numbers show a rather weak temperature 
dependence in the considered temperature range. The same trend is also seen at the 
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Figure 5. Partial static structure factors Sd/siq) (a, l3 = Si,0) for different 
densities at the temperature T = 2750 K. a) 5siSi{g), b) Sgioig), 'S'oo('?)- In 
panel d), Soo{q) for p = 4.2g/cm^ at three different temperatures is shown, 
indicating that the system has crystallized at T = 2750 K. 



two high densities, p = 3.7g/cm'^ and p = 3.9g/cm'^. Whereas the number of five-fold 
coordinated Si atoms keeps approximately constant around 50%, the coordination 
number 2310 = 6 further increases, whereas the tetrahedral coordination zsio — 4 
becomes less and less important with increasing density. At p = 3.9g/cm^, there are 
more oxygen atoms that are three-fold coordinated by Si atoms than bridging oxygens. 

The SiO coordination as a function of pressure is shown in Fig.[4]for T — 2100 K, 
2230 K and 3580 K. Although the qualitative behaviour of P(zsio) is similar at 
the different temperatures, there are significant quantitative differences between the 
coordination numbers at T = 3580 K and those at the two lower temperatures. The 
maximum of P{zsio = 5) shifts from about 30GPa at 3580 K to about 25GPa at 
the two lower temperatures. Moreover, at low temperatures, P{zsio — 4) decays 
more rapidly, accompanied by a faster increase of P(zsio — 6). Note that in 
recent simulation studies using both ab initio techniques and classical MD simulation 
[25l [26] , the pressure dependence of the SiO coordination has been only discussed 
for temperatures around 3500 K. However, Fig. |4] demonstrates that at least in the 
temperature range between 2000 to 4000 K the temperature dependence of P{zsio) 
has to be taken into account to allow for a quantitative comparison with experimental 
data, done at much lower temperatures. 
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Figure 6. Mean squared displacements for the silicon atoms for different densities 
at a) T = 2750 K and b) T = 2100 K. 



The inspection of the SiO and OSi coordination number distributions indicates 
that significant changes occur in the network structure for pressures above, say 10 GPa. 
This finding is similar to what has been observed experimentally for silica glasses. At 
pressures above about 20 GPa a network with five- and six-fold coordinated Si atoms is 
formed, similar to what has been recently found in other simulation studies [101 [^5 , 26J. 
It remains an open question whether in silica a first-order phase transition from a low 
density liquid to a high density liquid exists. Our data in Figs.[3]and[3]is also consistent 
with a gradual crossover from a low density phase (with zgio = 4) to a high density 
phase (with zgio — 6), occurring at low temperatures. 

To investigate the changes in the network structure on intermediate length scales, 
we now consider the partial static structure factors, defined as [18] 

^a/3(g) = ^^^(exp(iq.rM)) a,/3 = Si,0, (2) 
fe=i 1=1 
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with Na the number of atoms of type a, q the wavevector, and r^; — r; the 

distance vector between particle k and particle 

Fig. [5] shows the three Saf:i{q) at the temperature T = 2750 K for different 
densities. For the density 2.3g/cm^ the first peak [also called first sharp diffraction 
peak (FSDP)] is located at a wavenumber of about 1.6-1. 7 A~^. This peak is due 
to the tetrahedral order present in "low-density" silica. Its location is related to the 
length scale of two connected tetrahedra [30] . With increasing density, the first peak 
shifts to higher values of q and its amplitude decreases significantly. The behavior of 
the FSDP reflects the change of the structure from an open tetrahedral network at low 
densities to a more densely packed network at high densities, provided by SiOs and 
SiOg units as well as oxygen atoms that are three-fold coordinated by silicon atoms. 
Note that the partial structure factors Sapiq) are the required input for calculations 
in the framework of the mode-coupling theory of the glass transition [34], MCT, from 
which the dynamics of the melt can be predicted. A MCT calculation based on the 
Sai3 (q) of our simulation is presented elsewhere in this issue [35] . 

Figure [5jl shows Soo{q) for the density p = 4.2g/cm'^ at three different 
temperatures. At this high density, the system crystallizes at a temperature of 
T = 2750 K which is reflected by the occurrence of Bragg peaks in 5*00 (9)- At a 
temperature of T = 3250 K the system remains liquid; the structure factor in this case 
reminds one to that of a simple liquid, showing no sign of the FSDP feature. 

The structural changes in the network with increasing pressure are accompanied 
by an anomalous behavior of relaxation processes. This can be inferred from the 
mean squared displacements of a tagged particle of type a ~ Si, O, {rl.{t)) ~ 
(1/iVa) ^^^((rfc(t) — rfc(O))^). In Fig. [6| only (rgi(t)) is shown, since the mean 
squared displacement for oxygen exhibits a similar behaviour. At T = 2750 K, several 
time regimes can be clearly distinguished at the low density p = 2.21 g/ cm^, consisting 
of a ballistic regime at very short times {{r'^(t)) oc t'^), an intermediate regime with 
a maximum 0.2 ps, a plateau-like region and finally a linear regime, {r'^{t)) oc t^, 
indicating diffusive behaviour. Note that the maximum at 0.2 ps is both related to 
the so-called boson peak and to finite size effects present in the dynamics of silica 
(for details see Refs. [121 131] )• That this feature is not seen at high density might be 
related to the changes occurring in the network structure. From p — 2.21 g/cm^ to 
p = 3.5g/cm^, an acceleration of the diffusion dynamics can be inferred from Fig. [5^, 
followed by the opposite behaviour for a further decrease of the density. The same 
behaviour as for T = 2750 K is also seen for the lower temperature, T = 2100 K, 
although the effect of a maximal speed of diffusion is now more pronounced (Fig. [6)3) . 
One can also infer from Fig.[6l3 that the particles become more localized with increasing 
density, indicated by a decrease of the plateau's height with increasing density. 

The self-diffusion constants can be easily determined from the mean squared 
displacements using the Einstein relation. Da — lim(^oo(?'Q(i))/(6i). In Fig. [71 the 
self-diffusion constant for silicon, Dsi, along different isotherms is displayed as a 
function of density and pressure. Note that Dq, albeit slightly larger, exhibits the same 
behaviour. Along the different isotherms a maximum is found around about 3.5 g/cm"^, 
corresponding to a pressure of about 20 GPa (see Fig. \7]p) . Thereby, the diffusivity 
maximum becomes more pronounced with decreasing temperature. The behaviour of 
Dsi is closely related to that of the coordination numbers. Above about p = 20 GPa, 
the location of the diffusivity maximum, the number of five-fold coordinated Si atoms 
roughly shows a saturation, while the number of six-fold coordinated Si atoms further 
increases with increasing pressure. It seems that a SiO network with many five-fold 



MD simulation of amorphous silica under high pressure 



10 



10 

_ 10' 

f 10-' 
Q 10 

10 



■7 _ 



10' 



a) 4300K ^ 


.A ■■■ 


...A A--A.-A A 


A" 

* " 3580K 









2900K. 




-V — 


-'V V 


2750K 
















• 


a" , 




2100K * ^ 




•'' 






2580K .'^ 








2230K 






CHIK-SiOg 



2.0 



2.5 



3.0 

P (g/cm^) 



3.5 



4.0 



10' 
_ 10' 
I 10' 
Q 10 

10 



■7 



10' 



b) 



4300K 
3580K 



__2900K 

::::j»2750K 

2580K 
~^2230K 
*2100K 



CHIK-SiO. 



0.0 10.0 20.0 30.0 
P (GPa) 



40.0 



Figure 7. Sclf-difFusion constant of silicon as a function of density along different 
isotherms as a function of density, a), and as a function of pressure, b). 



coordinated Si atoms tends to have a faster diffusion dynamics tlian one with many six- 
fold coordinated Si atoms. Therefore, for pressures up to about 20 GPa the network 
dynamics is dominated by an increasing number of five-fold coordinated Si atoms, 
leading to an acceleration of diffusion. However, at pressures above about 20 GPa, the 
network dynamics is dominated by an increasing number of SiOe units, leading to a 
decrease of the diffusion constants with pressure. Thus, the diffusivity maximum is 
due to an interplay between five- and six-fold coordinated Si atoms. 

4. Conclusions 

Results of an extensive MD simulation have been presented to study structural and 
dynamic properties of liquid silica under (high) pressure. As an interaction potential 
a simple pair potential, the so-called CHIK potential, was used that leads to a reliable 
description of the equation of state. We have shown that the change of the network 
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structure with increasing pressure is related to an anomalous diffusion dynamics, 
characterized by a diffusivity maximum around 20 GPa. Elsewhere in this issue [3S] , 
a calculation in the framework of mode coupling theory based on the structural input 
of our simulation is presented that explicitly shows that the changes in the network 
structure are intimately related to the anomalous transport in amorphous silica. It 
is still an open question whether, at low temperatures, there is a first-order phase 
transition from a low-density to a high-density transition in silica. This question will 
be addressed in future studies with the CHIK potential. 
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